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INTRODUCTION 


In  order  to  enable  an  industrial  machine  with  primitive  computational  abil¬ 
ity  to  use  complicated  or  difficult  to  compute  functional  relationships 
repeatedly,  efficiently,  and  accurately,  it  is  necessary  to  supply  the  machine 
with  these  functional  relationships  as  sets  of  data  in  tabular  form.  It  is 
assumed  that  the  machine  can  deal  with  continuous,  piecewise  linear  functions 
(linear  splines).  A  graphics  tube  is  a  good  example.  Such  a  tube  can  draw  only 
straight  lines,  but  drawing  many  short,  connected  line  segments  can  represent  an 
arbitrary  curve  well.  In  order  to  represent  these  functions  most  accurately,  a 
nonuniform  mesh  must  be  used.  Finding  such  a  mesh  is,  in  principle,  a  very  dif¬ 
ficult  nonlinear  optimization  problem,  but  C.  de  Boor  (refs  1-3)  advocated  a 
general  method  by  which  the  mesh  can  be  found  quickly,  easily,  robustly  (and 
approximately)  without  any  recourse  to  optimization  methods!  We  present  herein 
a  robust  addition  to  de  Boor's  standard  method  which  improves  its  accuracy 
without  increasing  the  essential  complexity  of  his  algorithm. 

INTERPOLATORY  ERROR 

Let  i  be  the  linear  interpolant  of  function  f  on  a  subinterval  of  length  h. 
The  error  is  given  by 

f(t)  =  l(t)  +  e{t)  {M-^<t<u+^) 

Expand  e  in  a  Taylor  series  around  the  midpoint  (u)  of  the  subintsrval 

00  (i) 

e(t)  =  I  (t-u)^ 

i=0  ’  • 

Applying  the  two  boundary  conditions 

e(M  -  ^)  =  0  =  e(M  p 

ultimately  yields 


1 


^  u)  =  f  tilll  2i_ 

'2  (2i)!  ^2^  ^  ' 


00  (2-i-^l) 

/  (2i+l)!  '2'  '' 


Taking  the  first  two  terms  of  each  sum 


e("  +  M)  =  h»(t*-l) 

2  2“ 


2*0 


hat(t*-l) 


Letti ng 


one  has 


where 


and 


*  h*(t*-l) 


(5j 

2S.3.5  ^ 


{2+i) 

Pi  *  f(y)  /f”(M) 


>  U)  =  hMt*-l)(l  f  ht  +  -j2.  h*(t^O) 

+  hn(t"  +  n  +  0(h*)|  =  h*  (t*-l )  (l+S) 

S  =  a-iht  +  a2h*(t*+l)  +  a3h3t(t*  +  l)  +  0(h*) 


Pi  P2 

=  2T3  '  ®2  =  2*23 


P3 

®3  *  2»*325 


INTERPOLATORY  ERROR  NORM 

The  local  norm  of  the  error  on  a  subinterval  of  length  h  is  defined  by 

U+h/2 

"elln.h  =  (/  ,,,  le(t)  |  ^dtl^/" 

U-h/2 


2 


where  1  ^  n  <  oo  and  n  is  an  integer.  For  n  =  »,  we  have  the  maximum  error 


"elln.h  =  I  Li 


e(~  +  ji)  =  hMt*-l)(l+S) 

So  if  we  let  h  be  sufficiently  small  so  that  |  S  |  <  1  on  (-1,1),  we  have 


.  LnsiiV":) 


Since  only  the  even  terms  of  (1+S)'^  contribute  to  the  integral,  we  have 

llelln.h  =  (l-tM"EWUS)"dt 

where  Evd+Sli^  denotes  the  even  terms  of  (1+S)^. 


Hence, 


Ev(l+S)"  =  1  +  (J)a2hMl+t*)  +  (2)aih*t*  +  0(h*) 


Letting 


we  therefore  have 


^n,i  =  L 


/  (l-t*)^Ev(l+S)^dt 

0 


=  (l-tdf’(l+hMna2(l+td  +  aftM  +  0(h‘))dt 

.  In.o  nhMa2(In,o-^In,l)  *  "r  ailn.l)  *  0{h*) 


-  In,o(l+nhMa2(l  +  y'”)  *  "i"  ®1  7“"^  * 

*^n,o  ‘  ^n,o 

Using  integration-by-parts  on  and  solving  the  resulting  recursion  ult 


mately  yields 


2^"n! (2i)! (i+n) ! 

Ln,i  =  "~iTT2T+2n+nT'~ 
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from  which  we  conclude  that 


ahd 


Hence, 


-  fan^lTT" 

^n.l  -  ■■72n+3TT 


IdU  , 

In,o  2n+3 


(l-t*)^Ev(l+S)"dt 


and 


or 


where 


^®**n ,  h 


fafiJITT'  ai)  *  0(n<)) 


: . .  _ _ dh 


Hell, 


,,h  .  k  I  r'|„,  I  h!“/''|uhM2l5;|]  aj  .  5j5:3„  aj)  .  O(k-I) 


K  =  1 

2  '(2n+l)! ^ 

Using  Stirling's  approximation  to  the  factorial,  it  is  easy  to  show  that 


1 im  k  s  - 
n-oo  8 


Recalling  that 


we  finally  have 


Pi  P2 

ai  =  ^-5  and  32  =  2*~3 


Hell 


,2+1/n,,  ^  n;.2 


n-1  2 


n.h  =  M  f(u}  I  h  (1  *  24  <2nO  *  3(2F;;3T  ^ 


as  h-0,  where 


4 


and 


Pi 


1  f  (n!)*  i/n 

2  'T2n+iy!^ 


(2+i) 

=  f(4)  /f"{u) 


NORM  OF  ARBITRARY  FUNCTION 

The  local  norm  of  arbitrary  function  4>  over  a  subinterval  of  length  h  is 
defined  as 

M+h/2 

=  (/  U(t)  I  Pdt)^/P 

^  M-h/2 

where  p  >  0,  finite  and  real.  In  this  context,  we  allow  p  <  1  even  though 
Minkowski's  triangle  inequality  holds  only  for  p  ^  1. 

Expand  (()  in  a  Taylor  series  around  the  midpoint  of  the  subinterval 


where 


Now, 


but 


where 


and 


<(>(t)  =  X  ^7  (t-M)'' 

1=0  ''  • 


(i) 

Pi  *  (J>(M)/<()(M) 


+  4)  *  «<>{4)(1  +  S) 


S  =  I  ait’ 
i*l 


Pi  ,h  ’ 
=  il  ^2^ 


Hence,  letting  h  be  sufficiently  small  so  that  jS  |  <  1  on  (-1,1),  we  have 


(Id)  II 


p,h  =  ^  I  I  (l+S)Pdt  =  ^  I  tD(4)  1  Ev(l+S)Pdt 
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but 


hence 


We  therefore  have 


S  =  a^t  +  a2t*  +  a3t*  +  0(h*) 
£v(l+S)P  =  1  +  (j)a2t*  +  (2)81^* 


or 


as  h-0. 


Il<fttp,h  =  h  1  1  1  PtMa2  +  -5-  a^)  -  0(h-)dt 

=  h  1  1  P(1  +  IJ-  {P2  +  (p-l)Pi)  +  0(h-)) 


““'“p.h  “  I  «>(M)  I  (1  *  (P2  +  (P-l)Pi)  +  OCh-)) 


STANDARD  APPROXIMATION  TO 
Recalling  that 


and 


h'^'Plf'lp^h  •  I  f"(ul  I  (1  I-  5|(P2  *  (P-I)p{)  »  0(h*)) 


I  f..(m  I  ,i  .  ?i(5i25  pj  »  j.Sijj.  pj)  .  o(h-, 


24'2n+3 
.2+1/n 


we  multiply  the  first  equation  by  kh  and  subtract  from  the  second,  getting 


lleiln^h  = 


>  kh2*i/nif(p)  I  (^!(-  511-  P2  .  (I51|  -  p)p^  .  0(h*)) 


6n+9 


If  we  now  let  p  =  - — we  have 
2n+l 


^®*ln,h  *  *^**'^'''*n/ { 2n+l ) , h 

i-i.4+l/n  1  1  I  ,1  ,  n-*-!  ?(n*  +  14n+8 

^  kh  |f-(M)  I  (-J(-  --;5  P2  .  i2nll25nl9  ^1  ’  * 

=  kllfllp^h  ♦  kh^'"^^"’  !  f"(n)  !  (l-(-  ap2  1-  bPi)  i-  O(hM) 
=  kllf'’lln/,2n+l).h  * 


6 


For  n  =  1,2,  and  <»,  respectively,  we  have 


lleili.h  =  12  ‘l^"Wl/3,h  +  0(h») 

Ilell2,h  =  li^”»2/5,h  ^ 

2K30 

llelloo.h  =  §  »1=”Ml/2,h  +  Oth-*) 

STANDARD  ERROR  EQUIDISTRIBUTION  FOR  ANY  BANACH  NORM 

In  this  section,  we  justify  the  standard  method  of  error  equidistribution 
with  respect  to  any  Banach  norm.  The  global  norm  of  the  error  over  interval 
(a,D)  is 

llelln  =  {/  I  e(t)  I  "dtl^^" 


Hence ,  for  a  mesh  a  =  x-^  <  X2  <  ...  <  xn  =  b 

r,  N-1  Xi  +  i 


N-1 


I  /  ■  *  !  ei(.t)  Tdt  =  I  llelln, h, 
i  =  l  i=l 

Let  single  bars  around  the  error  denote  the  standard  aoproximation  to  the 
error  norm  and  analogously  define 


but 


where 


Hence,  letting 


we  have 


e 


n 

n 


N-1 

li 


j=l 


e 


n 

n,hj 


><jt-l 

n,h.  =  kHf'il  n.  =  •‘(Z  I  I  '"dti'/P 

J  P'"l  X, 


P 


n 

2n+i 


Xj+i 

ip,hi  =  /  !f"(t: 

J  Xj 


"dt 


n  N-1 

e  I  n  =  •«"  I 

j=l 


2n+l 
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We  win  refer  to  the  integrals  the  standard  or  de  Boor  integrals. 

It  follows  trivially,  using  Leibnitz's  rule,  that 


implies  that 


a  .  if' 

I  e  I  =0  1  <  i  <  N 

°x.j  n 


Ip,h.i_l  ■  ^P'^i  1  <  ■>  <  N 

Hence,  the  condition  Ip  f,  =  constant  determines  the  mesh  which  minimizes  the 
standard  global  approximation  to  (lelln- 

For  a  linear  spline  approximation  to  f",  it  is  a  fairly  simple  (see 
COMPUTATION)  matter  to  find  the  mesh  for  which  the  de  Boor  integrals  are 
constant . 


CONVERGENCE  OF  STANDARD  METHOD 

Recall  that 

llelln,h  =  '«ll^""p,h  ^  54  f’'()J)  1  (-ap2+bo|)  + 


Letting 


F  =  -aflylflMl  +  bfi 
we  have  the  following  one  term  approximation  to  the  difference  between  Helln  ^ 
and  !  el 

4+1/n 


but 


Therefore,  we  also  have 


"®"n,h  -  !  ®  1  n,h  = 


lleHn^h  =  f"(u) 


''®''n,h  ■  '  ®  I  n,h  Fh* 

iieH";;  “  24rf"(4)T* 


but  also 

nf'll  ^  =1  f"(p)  1 

P ,  h 
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hence , 

h  a  (, - --V)  =  (i - = - - 

T^”(m)T’  T^'‘{m)T’  !f>'(4)|P 

In  addition,  for  the  correct  mesh 

Ip,h  =  N-i 

hence , 

h  ,  - - 

(M-l)l  f(M)  I  P 

and  therefore, 

Hell  L^'lsl  K  FT* 

n^h _ _ _ __B _ 

"«“n,h  24l  f'(u)  1  2 

This  tells  us  that  the  relative  difference  between  llellp,^^  and  |  e  |  ,-|^h 

as  N-^,  which  means  that  tne  standard  method  works  better  and  better 

iHelln.h  P®  more  nearly  constant)  as  N  gets  larger  and  larger.  This  is  all 

true  however,  with  the  proviso  that 

_  __  F _ 

!  f(p)  I 

IS  bounded  throughout  the  region  of  interest.  It  stands  to  reason,  therefore, 
that  the  standard  method  will  perform  worst  where  f"  is  not  bounded  away  from 
zero . 

IMPROVED  APPROXIMATION  TO  leln.h 
Recall  that 

M  '"(u)  I  .  51(2225  "Z  *  3nR;5)  “i>  * 

and 

llf"llq,h  =  !  (1  ^5(P2  +  (q-l)Pi)  *  0(h*)) 

Multiplying  h  by  r  in  the  second  equation,  we  have 

r'^/'^tlf-llq  rh  =  "  rMq-l)Pi)  +  0(h*)) 


9 


Multiplying  this  equation  by  khQ  gives  us 


f”(U)  I  (1  +  +  rMq-l)Pi)  +  0(h*; 


q,  rh 


Now,  in  order  to  make  this  equation  look  as  much  like  the  very  first  one  as 
possible,  we  set 

n-l 


r*  =  Ql?-  ,  p2{q_l)  5. 

2n+3  '  '  3(2n+3: 


and 


Q  4.  1  =  2  - 

q  n 


Solving  for  r,  q,  and  Q,  we  have 


,n+2  1/2 

''  “  -2043 


4n+§ 

3n+6 


and 


Q  = 


5n*+8n+5 

4n*+5n 


A  simple  subtraction  then  gives  us  an  improved  approximation  to  llell^^  ^ 

II  ei' 


n ,  h 


where  before,  we  had 

'"'n.h  =  ''''"'p.h  •  =  I  '  I  „,h  •  OO'*’''") 

It  must  be  mentioned  however,  that  although  this  improved  approximation  is 
asymptotically  more  efficient,  no  such  approximation  can  be  uniformly  superior 
in  all  cases.  Bearing  this  in  mind,  we  dispense  with  approximations  on  all 
subintervals  not  having  f"  bounded  away  from  zero  and  instead  use  the  exact 
error 

^  ^  x—x  ■  ^i  +  1  ^ 

®i(^)  =  /  /  f"(u)dudt  -  - - —  /  /  f"{u)dudt 

^1  ^i  ^i+l~>^i  ^i  ^i 
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COMPUTATION 

In  actual  computation,  we  assume  the  existence  of  a  piecewise  linear 
approximation  to  |  f"  j  .  The  mesh  over  which  this  function  is  defined  is 
referred  to  as  the  "original"  mesh.  In  order  to  deal  with  the  standard  and 
improved  asymptotic  integral  approximations  to  the  local  error  norm,  we  will 
need  to  deal  with  integrals  of  the  form 

c+1 

L  =  /  X(t)i"/f'dt 
c 

where  A  is  a  nonnegative  linear  function  with  slope  s 

A(t)  =  A(c)  +  s(t-c) 

wi  th 


A(t)>0  for  c<t<c+i 
and  where  m  and  n  are  arbitrary  positive  integers. 

In  the  following,  let 

a  =  X(c)^^'^ 
and 

k  .  .  „k+l  „k+l 

Sk  =  I  >>  s'"’  = 

Pirst.  we  need  to  compute  l.  as  a  function  of  £ 

,  _  ■•f’Sm+n-i 

TmInTc"”" 

(m+n)Sn_^ 

where 

0  =  (X(c)  4- 

Second,  we  need  to  compute  I  as  a  function  of  L 

L(m+n)Sn_i 

,  =  - D_l  =  B(L) 

^®m+n-l 

where 

0  =  (Xlc)"*/"*^  t  (?  .  l)sL)l/<'"^"’ 
n 
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A  and  B  are  therefore  inverse  functions,  i.e., 

A(B(x))  =  X  =  B{A(x)) 
or 


A-’  B  and  B'’  =  A 

Now  let  values  of  u  denote  the  original  mesh  and  let  g  be  the  piecewise  linear 
interpolant  to  the  (u^.l  f^"  ]  )  data. 

Define  the  integral 

G(x)  =  /  g(t)'"/"dt 

‘^1 


Now  if  u-i  $  X  $  u-j+i, 

Ui  ,  Ui+X-U,  , 

G{x)  =  /  g(t)'"/"dt  /  g^itl^^^dt 

Ui  Ut 


=  G(Ui)  +  L 

where  X  =  g-j ,  c  =  u^ ,  and  i  =  x-u^.  Hence, 

G(x)  =  G(ui )  A(x-Ui ) 

explicitly  defines  G  for  all  x  in  the  domain  of  interest. 

In  order  to  get  the  standard  mesh,  we  will  also  have  to  compute  the  inverse 
of  G  (only  for  m/n  =  p) . 

B(G(x)  -  G(Ui))  =  8(A(x-Ui))  =  x-u. 

Hence , 


X  =  Ui  +  8(G(x)  -  G(Ui ) ) 
but  if  G{x)  =  y,  then  x  =  G'Mr).  Therefore, 

G-’(r)  =  Ui  +  B(y  -  G(ui)) 


for 


and  provided 


G(Ui )  <  r  <  G(Ui+i ) 


G(Ui)  *  G(Ui+i) 
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Define 


li  =  G(Xt+i)  -  G(xi)  =  /  g(t)Pdt 

^i 

ere  x  is  the  standard  or  improved  mesh,  obtained  by  prescribing  values  for  the 
I's.  The  standard  method  prescribes 

G(Xn) 

li  =  const  =  (1  <  i  ^  N) 

For  the  improved  mesh,  the  I's  will  vary,  but  the  mesh  is  still  obtained  in  the 
standard  way.  Since 

G(x.i  4.^  )  =  G(x.j  )  +  I  i 

we  have  immediately  that 

Xi  +  i  =  G-'(G(Xi)  +  I.i)  i  =  1,2 . N-2 


ALGORITHM 

Let  *  denote  a  standard  or  improved  mesh  and  **  denote  the  succeeding 
improved  mesh.  We  have  seen  that  the  main  contributor  to  the  ratios 
llelln,h**/llelln,h*  and  i  e  |  n,h’^/1  e  I  n,h*  is 

,h!!,2+l/n 
h*  f’'(M*) 

We  therefore  have  the  approximate  asymptotic  relation 


I  e  I  r,.h*  “  lle«n,h* 


But  we  would  like  lie II p  to  be  constant,  hence  we  have  the  proportionality 


n,  h 


n,h’ 


'>e''n,h* 


or 


*  P 


T  **  ®  '  fi-fi  . 
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We  calculate  the  I's  accordingly  and  multiply  them  by  the  appropriate  constant 


to  get 


N-1 

I  ip.n* 

i  =  l  ■' 


'n-I" 


The  quantities  llellp^h*  are  computed  either  from  the  improved  asymptotic 
approximation  or  exactly  (relative  to  the  original  data)  depending  on  whether  or 
not  f"  is  bounded  away  from  zero  on  the  subinterval  in  question.  It  is  impor¬ 
tant  to  note  that  this  approximate  relation  between  the  *  and  **  meshes  can  lead 
to  exact  convergence  (rapidly)  to  the  minimax  mesh.  If  the  *  mesh  ^  the 
mimmax  mesh  (llelln,h*  =  constant),  then  the  de  Boor  integrals  (Ip,h) 
mesh  will  be  no  different  from  those  on  the  *  mesh. 

The  practical  convergence  properties  of  this  algorithm  are  as  follows.  If 
f"  is  well  bounded  away  from  zero,  the  standard  de  Boor  method  gives  impeccable 
results  without  any  iteration.  If  f”  is  not  bounded  away  from  zero,  convergence 
to  a  virtually  perfect  minimax  mesh  can  easily  occur  in  only  two  iterations.  A 
few  iterations  may  be  needed  in  the  presence  of  multiple  inflection  ooints. 

In  any  case,  even  the  very  first  iteration  improves  the  mesh  markedly. 
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